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Abstract. A widely used electrostatics model in the biomolecular modeling com- 
munity, the nonlinear Poisson-Boltzmann equation, along with its finite element ap- 
proximation, are analyzed in this paper. A regularized Poisson-Boltzmann equation is 
introduced as an auxiliary problem, making it possible to study the original nonlinear 
equation with delta distribution sources. A priori error estimates for the finite element 
approximation are obtained for the regularized Poisson-Boltzmann equation based on 
certain quasi-uniform grids in two and three dimensions. Adaptive finite element ap- 
proximation through local refinement driven by an a posteriori error estimate is shown 
to converge. The Poisson-Boltzmann equation does not appear to have been previously 
studied in detail theoretically, and it is hoped that this paper will help provide molecular 
modelers with a better foundation for their analytical and computational work with the 
Poisson-Boltzmann equation. Note that this article apparently gives the first rigorous 
convergence result for a numerical discretization technique for the nonlinear Poisson- 
Boltzmann equation with delta distribution sources, and it also introduces the first prov- 
ably convergent adaptive method for the equation. This last result is currently one of 
only a handful of existing convergence results of this type for nonlinear problems. 
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1. Introduction 

In this paper, we shall design and analyze finite element approximations of a widely 
used electrostatics model in the biomolecular modeling community, the nonlinear Poisson- 
Boltzmann equation (PBE): 

- V ■ (eVn) + sinh(M) = ^ qi6i in M'^, ci = 2, 3, (1.1) 

i=l 

where the dielectric e and the modified Debye-Hiickel parameter R are piecewise con- 
stants in domains Vim (the domain for the biomolecule of interest) and ^Is (the domain for 
a solvent surrounding the biomolecule), and Si := S{x—Xi) is a Dirac distribution at point 
Xi. The importance of (1.1) in biomolecular modeling is well-established; cf. [14, 43] 
for thorough discussions. Some analytical solutions are known, but only for unrealistic 
structure geometries, and usually only for linearizations of the equation; cf. [29] for a 
collection of these solutions and for references to the large amount of literature on ana- 
lytical solutions to the PBE and similar equations. The current technological advances 
are more demanding and require the solution of highly nonlinear problems in compli- 
cated geometries. To this end, numerical methods, including the finite element method, 
are widely used to solve the nonlinear PBE [29, 30, 5, 6, 44, 19, 56]. 

The main difficulties for the rigorous analysis and provably good numerical approx- 
imation of solutions to the nonlinear Poisson-Boltzmann equation include: (1) Dirac 
distribution sources, (2) exponential rapid nonlinearities, and (3) discontinuous coeffi- 
cients. We shall address these difficulties in this paper. To deal with the 5 distribution 
sources, we decompose u as an unknown function in and a known singular function, 
namely, 

u = u + G, with G = ^ Gi, 

i=l 

where Gi is the fundamental solution of —Sm^Gi = qi5i in W'-. Substituting this de- 
composition into the PBE, we then obtain the so-called regularized Poisson-Boltzmann 
equation (RPBE): 

-V ■ {eVu) + sinh(u + G) = V ■ ((e - £™)VG) in = 2, 3. 

The singularities of the 5 distributions are transferred to G, which then exhibits degen- 
erate behavior at each {xi} C Vtm- At those points, both sinh G{xi) and VG{xi) exhibit 
blowup. However, since G is known analytically, one avoids having to build numerical 
approximations to G. Moreover, both of the coefficients k and e — Em^^ zero inside Vim 
where the blowup behavior arises. Due to this cutoff nature of coefficients, we obtain a 
well-defined nonlinear second-order elliptic equation for the regularized solution u with 
a source term in . We will show that it also admits a unique solution u E H^, even 
though the original solution u ^ due to the singularities present in G. 

Singular function expansions are a common technique in applied and computational 
mathematics for this type of singularity; this type of expansion has been previously pro- 
posed for the Poisson-Boltzmann equation in [58] and was shown (empirically) to allow 
for more accurate finite difference approximations. In their work, the motivation for the 
technique was the poor discrete approximation of arbitrarily placed delta distributions 
using only the fixed corners of uniform finite difference meshes. In the present work, our 
interest is in developing finite element methods using completely unstructured meshes, so 
we are able to place the delta distributions precisely where they should be and do not have 
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this problem with approximate delta function placement. Our motivation here for consid- 
ering a singular function expansion is rather that the solution to the Poisson-Boltzmann 
equation is simply not smooth enough to either analyze or approximate using standard 
methods without using some sort of two-scale or multiscale expansion that represents the 
nonsmooth part of the solution analytically. In fact, it will turn out that expanding the 
solution into the sum of three functions, namely, a known singular function, an unknown 
solution to a linear auxiliary problem, and an unknown solution to a second nonlinear 
auxiliary problem, is the key to establishing some fundamental results and estimates for 
the continuous problem and is also the key to developing a complete approximation the- 
ory for the discrete problem as well as provably convergent nonadaptive and adaptive 
numerical methods. 

Starting with some basic results on existence, uniqueness, and a priori estimates for the 
continuous problem, we analyze the finite element discretization and derive discrete ana- 
logues of the continuous results to show that discretization leads to a well-posed discrete 
problem. Using maximum principles for the continuous and discrete problems, we derive 
a priori L°°-estimates for the continuous and discrete solutions to control the nonlinear- 
ity, allowing us to obtain a priori error estimates for our finite element approximation of 
the form 

\W - Uh\\i < inf \\u - Vh\\i, 

where is the linear finite element subspace defined over quasi-uniform triangulations 
with a certain boundary condition, and Uh is the finite element approximation of u in V^. 
The result is quasi-optimal in the sense it implies that the finite element approximation 
to the RPBE is within a constant of being the best approximation from the subspace V^. 
After establishing these results for finite element approximations, we describe an adap- 
tive approximation algorithm that uses mesh adaptation through local refinement driven 
by a posteriori error estimates. The adaptive algorithm can be viewed as a mechanism for 
dealing with the primary remaining difficulty in the RPBE, namely, the discontinuities 
of the coefficients across the interface between the solvent and the molecular regions. 
Finally, we shall prove that our adaptive finite element method will produce a sequence 
of approximations that converges to the solution of the continuous nonlinear PBE. This 
last result is one of only a handful of existing results of this type for nonlinear elliptic 
equations (the others being [24, 48, 15]). 

The outline of this paper is as follows. In section 2, we give a brief derivation and 
overview of the Poisson-Boltzmann equation. In section 3, we derive a regularized form 
of the Poisson-Boltzmann equation by using a singular function expansion. In section 4, 
we give some basic existence and uniqueness results for the RPBE. In section 5, we de- 
rive an a priori L°° -estimate for the continuous problem. After introducing finite element 
methods for the RPBE, in section 6 we derive an analogous a priori L°°-estimate for 
the discrete problem, and based on this we obtain a quasi-optimal a priori error estimate 
for the finite element approximation. In section 7, we describe the adaptive algorithm, 
present an a posteriori error estimate, and prove a general convergence result for the al- 
gorithm. In the last section, we summarize our work and give further remarks on the 
practical aspects using results in the present paper. 

2. The POISSON-BOLTZMANN EQUATION 

In this section we shall give a brief introduction to the nonlinear Poisson-Boltzmann 
equation. A detailed derivation can be found in [47, 29]. 
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The nonlinear PBE, a second-order nonlinear partial differential equation, is funda- 
mental to Debye-Hiickel continuum electrostatic theory [22]. It determines a dimen- 
sionless potential around a charged biological structure immersed in a salt solution. The 
PBE arises from the Gauss law, represented mathematically by the Poisson equation, 
which relates the electrostatic potential $ in a dielectric to the charge density p: 

-V ■ (£V<I>) = p, 

where e is the dielectric constant of the medium and here is typically piecewise constant. 
Usually it jumps by one or two orders of magnitude at the interface between the charged 
structure (a biological molecular or membrane) and the solvent (a salt solution). The 
charge density p consist of two components: p = Pmacro + Pion- For the macromolecule, 
the charge density is a summation of 5 distributions at A^^ point charges in the point 
charge behavior, i.e.. 



Pmacro('^) ^ ^ Qi^j^-^ -^i): Qi 



where > is the Boltzmann constant, T is the temperature, Cc is the unit of charge, 
and Zi is the amount of charge. 

For the mobile ions in the solvent, the charge density pion cannot be given in a de- 
terministic way. Instead it will be given by the Boltzmann distribution. If the solvent 
contains N types of ions, of valence Zi and of bulk concentration Cj, then a Boltzmann 
assumption about the equilibrium distribution of the ions leads to 



N 

Pi 



For a symmetric 1 : 1 electrolyte, N = 2, Ci = cq, and Zi = (—1)*, which yields 

Pion = -2coec sinh 



(—)■ 



We can now write the PBE for modeling the electrostatic potential of a solvated bi- 
ological structure. Let us denote the molecule region by ilm C R.'^ and consider the 
solvent region = K'^V^m- We use u to denote the dimensionless potential and to 
denote the modified Debye-Hiickel parameter (which is a function of the ionic strength 
of the solvent). The nonlinear Poisson-Boltzmann equation is then 



-V ■ {eVu) + sinh(M) = ^ qiSi in R'^, (2. 1) 



i=l 



u{oo) = 0, (2.2) 



Em if X E ^Im, and ^ - if x G ^1^, 



where 



Es if X E ^Is, 1 y/^sf^ > if X E ^Is- 

It has been determined empirically that ~ 2 and Es ~ 80. The structure itself (e.g., 
a biological molecule or a membrane) is represented implicitly by e and R, as well as 
explicitly by the A^^, point charges qi = ZiCc at the positions Xj. The charge positions 
are located in the strict interior of the molecular region fi^. A physically reasonable 
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mathematical assumption is that all charge locations obey the following lower bound on 
their distance to the solvent region for some a > 0: 

k — > (T Vx G Qs, i = I, ■ ■ ■ , N^- (2.3) 

In some models employing the PBE, there is a third region Vti (the Stem layer [11]), a 
layer between Vtm and VLg. In the presence of a Stern layer, the parameter o in (2.3) 
increases in value. Our analysis and results can be easily generalized to this case as well. 

Some analytical solutions of the nonlinear PBE are known, but only for unrealistic 
structure geometries and usually only for linearizations of the equation; cf. [29] for a 
collection of these solutions and for references to the large amount of literature on an- 
alytical solutions to the PBE and similar equations. However, the problem is highly 
nonlinear. Surface potentials of the linear and the nonlinear PBE differ by over an order 
of magnitude [44]. Hence, using the nonlinear version of the PBE model is fundamen- 
tally important to accurately describe physical effects, and access to reliable and accurate 
numerical approximation techniques for the nonlinear PBE is critically important in this 
research area. 

We finish this section by making some remarks about an alternative equivalent formu- 
lation of the PBE. It is well known (cf. [47, 29]) that the PBE is formally equivalent to 
a coupling of two equations for the electrostatic potential in different regions Vt^ and Vtg 
through the boundary interface. This equivalence can be rigorously justified. Inside Vtrn, 
there are no ions. Thus the equation is simply the Poisson equation 

-V ■ (emVw) = ^ qi5i in fi^. 

i=l 

In the solvent region Vtg, there are no atoms. Thus the density is given purely by the 
Boltzmann distribution 

—V ■ (^s Vm) + sinh(-u) = in Vtg- 

These two equations are coupled together through the boundary conditions on the inter- 
face r := dVljn = dVLg n Vim'- 

[u\y = 0, and 

where [/]|r = Yim.t^Q f{x + tnr) — f{x — tnr), with nr being the unit outward normal 
direction of interface T. We will assume T to be sufficiently smooth, say, of class C^. 

Solving the individual subdomain systems and coupling them through the boundary, 
in the spirit of a nonoverlapping domain decomposition method, is nontrivial due to the 
complicated boundary conditions and subdomain shapes. Approaches such as mortar- 
based finite element methods to solve the coupled equations for linear or nonlinear PBE 
can be found in [19, 51]. 

3. REGULARIZATION of the CONTINUOUS PROBLEM 

In this section, we shall introduce a regularized version of the nonlinear PBE for both 
analysis and discretization purposes. We first transfer the original equation posed on 
the whole space to a truncated domain using an artificial boundary condition taken from 
an approximate analytical solution. Then we use the fundamental solution in the whole 
space to get rid of the singularities caused by 6 distributions. We shall mainly focus on 
more difficult problems in three dimensions. Formulation and results in two dimensions 
are similar and relatively easy. 



du 
drir 



6 



L. CHEN, M. HOLST, AND J. XU 



Let i7 C with a convex and Lipschitz-continuous boundary d^l, and C fi. 
In the numerical simulation, for simplicity, we usually choose to be a ball or cube 
containing a molecule region. The solvent region is chosen as f2s fl f2 and will be still 
denoted by On we choose the boundary condition u = g, with 
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^2 \ ^ ^-K\x-Xi\ 

kBTj ^ es\x - Xi\ ' 
' 1=1 



The boundary condition is usually taken to be induced by a known analytical solution 
to one of several possible simplifications of the linearized PBE. Far from the molecule, 
such analytical solutions provide a highly accurate boundary condition approximation 
for the general nonlinear PBE on a truncation of M^. For example, (3.1) arises from the 
use of the Green's function for the Helmholtz operator arising from linearizations of the 
Poisson-Boltzmann operator, where a single constant global dielectric value of Eg is used 
to generate the approximate boundary condition. (This is the case of a rod-like molecule 
approximation; cf. [29].) Another approach to handling the boundary condition more 
accurately is to solve the PBE with boundary conditions such as (3.1) on a large (with 
a coarse mesh) and then solve it in a smaller Q (with a fine mesh) with the boundary 
condition provided by the earlier coarse mesh solution. The theoretical justification of 
this approach can be found at [28] using the two-grid theory [53]. We are not going to 
discuss more on the choice of the boundary condition in this paper. 
Employing (3.1) we obtain the nonlinear PBE on a truncated domain: 

-V ■ (eVu) + sinh(M) = ^ q,Si in Q, (3.2) 

i=l 

u = g on dVL. (3.3) 

This is, in most respects, a standard boundary-value problem for a nonlinear second-order 
elliptic partial differential equation. However, the right side contains a linear combina- 
tion of 5 distributions, which individually and together are not in H^^iVt); thus we cannot 
apply standard techniques such as classical potential theory. This has at times been the 
source of some confusion in the molecular modeling community, especially with respect 
to the design of convergent numerical methods. More precisely, we will see shortly that 
the solution to the nonlinear Poisson-Boltzmann equation is simply not globally smooth 
enough to expect standard numerical methods (currently used by most PBE simulators) 
to produce approximations that converge to the solution to the PBE in the limit of mesh 
refinement. 

In order to gain a better understanding of the properties of solutions to the nonlinear 
PBE, primarily so that we can design new provably convergent numerical methods, we 
shall propose a decomposition of the solution to separate out the singularity caused by 
the 5 distributions. This decomposition will turn out to be the key idea that will allow 
us to design discretization techniques for the nonlinear PBE which have provably good 
approximation properties and, based on this, also design a new type of adaptive algorithm 
which is provably convergent for the nonlinear PBE. 

We now give this decomposition. It is well known that the function 

solves the equation 

-V ■ {smVGi) = q^6i in Ml 
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We thus decompose the unknown u as an unknown smooth function u and a known 
singular function G: 

u = u + G, 

with 

G = J2G^■ (3.4) 

i=l 

Substituting the decomposition into (3.2), we then obtain 

-V ■ (eVu) + sinh(M + G) = V • ((e - em)VG) in n, (3.5) 

u = g — G on dVL, (3.6) 

and call it the RPBE. The singularities of the 5 distribtuions are transferred to G, which 
then exhibits degenerate behavior at each {xj} C Vtm- At those points, both sinh G'(xj) 
and VG(xj) exhibit blowup. However, since G is known analytically, one avoids having 
to build numerical approximations to G. Moreover, both of the coefficients R and e — 
Em are zero inside Vtm^ where the blowup behavior arises. Due to this cutoff nature of 
coefficients, the RPBE is a mathematically well defined nonlinear second-order elliptic 
equation for the regularized solution u with the source term in H^^. We give a fairly 
standard argument in the next section to show that it also admits a unique solution u G 
H^, even though the original solution u ^ due to the singularities present in G. In 
the remainder of the paper we shift our focus to establishing additional estimates and 
developing an approximation theory to guide the design of convergent methods, both 
nonadaptive and adaptive. 

Before moving on, it is useful to note that, away from {xi], the function G is smooth. 
In particular, we shall make use of the fact that G G C°^(l],) n C~(r) n G°^{dn) 
in the later analysis. Also, a key technical tool will be a further decomposition of the 
regularized solution u into linear and nonlinear parts, u = v} + m", where v! satisfies 

-V ■ (eVu') = V ■ ((£ - e^)VG) in Vt, (3.7) 

= on dVL, (3.8) 

and where satisfies 

-V ■ (eVu") + R^ sinh(u" + + G) = in O, (3.9) 

yJ' = g-G ondVt. (3.10) 

4. EXISTENCE AND UNIQUENESS 

In this section we shall discuss the existence and uniqueness of the solution of the 
continuous RPBE. The arguments we use in this section appear essentially in [29], ex- 
cept there the PBE was artificially regularized by replacing the delta distributions with 
i/^^-approximations directly rather than being regularized through a singular function 
expansion. 

We first write out the weak formulation. Since AG = away from {xi], through 
integration by parts we get the weak formulation of RPBE: Find 

ueM ■={v e H\n) I e^ e"^ G L'^i^s), and v = g - G on dn} 

such that 

A{u,v) + {B{u),v) + {fG,v) = yveH',{n), (4.1) 

where 

• A{u,v) = (eVw, Vm), 
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• {B{u),v) = {k^ smh(n + G),v), and 

• {fG,v) = J^{e-ejVG-Vv. 
Let us define the energy on M: 

E{w) = / ^|Vwp + K^cosh(iy + G) + (/g,w). 

It is easy to cliaracterize tlie solution of (4.1) as the minimizer of the energy. 
Lemma 4.1. Ifu is the solution of the optimization problem, i.e., 

E{u) = inf E{w), 

then u is the solution o/(4. 1). 

Proof. For any v E Hq{Q) and any t G M, the function F(t) = E{u + tv) attains the 
minimal point at t = 0, and thus -F'(O) = 0, which gives the desired result. □ 

We now recall some standard variational analysis on the existences of the minimizer. 
In what follows we suppose 5 is a set in some Banach space V with norm || ■ ||, and J{u) 
is a functional defined on S. S is called weakly sequential compact if, for any sequence 
{uk} C S, there exists a subsequence {ufc^} such that Uk^ ^ u E S, where ^ stands for 
the convergence in the weak topology. For any u, if J{uk) — )■ J{u), we say J is 

weakly continuous at m; if 

J{u) < liminf J{uk), 

we say J is weakly lower semicontinuous (w.l.s.c.) at u. The following theorem can be 
proved by the definition easily. 

Theorem 4.2. // 

1 . S is weakly sequential compact, and 

2. J is weakly lower semicontinuous on S, 
then there exists u E S such that 

J{u) = inf J{w). 

We shall give conditions for the weakly sequential compactness and weakly lower 
semicontinuity. First we use the fact that a bounded set in a reflexive Banach space is 
weakly sequential compact. 

Lemma 4.3. One has the following results: 

1. The closed unit ball in a reflexive Banach space V is weakly sequential compact. 

2. 7/'lim||„||_j.oo J{'v) = oo, then 

inf J(w) = inf J(w). 

The next lemma concerns when the functional is w.l.s.c. The proof can be found 
at [57]. 

Lemma 4.4. If J is a convex functional on a convex set S and J is Gateaux differentiable, 
then J is w.l.s.c. on S. 

Now we are in the position to establish the existence and uniqueness of solutions to 
the RPBE. 

Theorem 4.5. There exists a unique u E M G H^{Q) such that 

Eiu) = inf E(w). 

w&M 
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{fG,v)<e,\\VG\\nA^v\\n^ < || VG||^ + v^||2 



Proof. It is easy to see E{w) is differentiable in M with 

{DE{u),v) = A{u,v) + {B{u),v) + {f^v). 

To prove the existence of the minimizer, we need only to verify that 

(1) M is a convex set, 

(2) E is convex on M, and 

(3) Um||„||^^oo-E(f) = oo. 

The verification of (1) is easy and thus skipped here. (2) comes from the convexity 
of functions x"^ and cosh(x). Indeed E is strictly convex. (3) is a consequence of the 
inequality 

Eiv)>Ci6,R)\\v\\l + CiG,g), (4.2) 
which can be proved as following. First, by Young's inequality we have for any S > 

1 

6' 

Since cosh(x) > 0, we have then > K)||Vt;|p — (1/(5)||VG||^^, where we can 
ensure C{e, R) > if 5 is chosen to be sufficiently small. Then using norm equivalence 
on M, we get (4.2). The uniqueness of the minimizer comes from the strict convexity of 
E. □ 

5. Continuous a priori L°°-estimates 

In this section, we shall derive a priori -estimates of the solution of the RPBE. The 
main result of this section is the following theorem. 

Theorem 5.1. Let u be the weak solution of RPBE in H^{Q). Then u is also in L°°{Q). 

Note that we cannot apply the analysis of [31, 32] directly to the RPBE, since the 
right side fc G H'^(^l) and does not lie in L°°(fi) as required for use of these results. 
We shall overcome this difficulty through further decomposition of u into linear and 
nonlinear parts. 

Let u = + m", where E Hq^VI) satisfies the linear elliptic equation (the weak 
form of (3.7)-(3.8)) 

Aiu',v) + {fG,v) = WveH'.in) (5.1) 

and where m" G M satisfies the nonlinear elliptic equation (the weak form of (3.9)- 
(3.10)) 

A{u'',v) + {B{u'' + u^),v) = ^veH^iQ). (5.2) 

Theorem 5.1 then follows from the estimates of and in Lemmas 5.2 and 5.3; cf. 
(5.3) and (5.4). 

Lemma 5.2. Let be the weak solution o/(5.1). Then 

G L°°(fi). (5.3) 

Proof. Since AG = in Vtg^ using integral by parts we can rewrite the functional fc as 



5 

r 



where [e] = Eg — Sm is the jump of e at the interface. We shall still use fc to denote the 
smooth function [e]^ on F. 
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It is easy to see that the linear equation (5.1) is the weak formulation of the elhptic 
interface problem 

du'' 



-V ■ (eVw') = in = 0, 



dn 



fc on r, and m = on dVl. 



Since fc G C°°(T) and T G C^, by the regularity result of the elliptic interface problem 
[4, 12, 20,41], we have G H'^{nrn)^H'^{^s)^H^{^)- In particular by the embedding 
theorem we conclude that v} E L^{Vl). □ 

To derive a similar estimate for the nonlinear part m", we define 

a' = argmax ( sinh(c + sup {v} + G)) < ) , a = min ( a', inf — G) ) , 

/3' = argmin ( sinh(c + inf (u' + G)) > ) , /3 = max ( sup(5f - G) ) . 
The next lemma gives the a priori -estimate of m". 

Lemma 5.3. Let Z^e weak solution of (5.2). Then a < u" < (3, and thus 

G (5.4) 

Proof. We use a cutoff-function argument similar to that used in [3 1 ] . Since the boundary 
condition g — G E G°^{d^l), we can find au^ E H^(^l) such that ud = g — G on dVl in 
the trace sense, or more precisely 

Tud = g-G, 

where T : i7 i— )■ dVl is the trace operator. Then the solution can be written u"^ = + uq, 
with no G H^i^). Let = - /3)+ = max(u" - /3, 0) and = («" - a)" = 
min(M" — a, 0). Then from 

< 0= K -/?)+ = (md + Mo -/?)+< (no -/?)+ + 
> = (m" — a)^ = {ud + uq — a)^ > {ud — a)" + , 

and 

< T0 < T(nB - + T< = 0, 
> T0 > T(ud - a)" +Tno = 0, 

we conclude that both 0, G ifo(i7). Thus for either = or = 0, we have 

(£Vm", V0) + (k^ sinh(M" + + G), 0) = 0. 

Note that > in and its support is the set 3^ = {x G | > (3}. On 3^, we have 

sinh(M" + + G) > sinh f + inf (m' + G)) > 0. 
Similarly, < in i7 with support %Q\.y_= {x E Vl\u^{x) < a]. On X we now have 
sinh(M" + + G) < sinh (a' + inf (m' + G)) < 0. 

Together this implies 

> (eVw", V0) = (eVK - V0) = ell V0f > 
for either = or = 0. Using the Poincare inequality we have finally 

O<||0||<||V0||<O, 
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giving = 0, again for either = or = 0. Thus a < < /3 in f2. □ 

6. Finite element methods for the regularized Poisson-Boltzmann 

equation 

In this section we shall discuss the finite element discretization of RPBE using linear 
finite element spaces and prove the existence and uniqueness of the finite element 
approximation u^. Furthermore, under some assumptions on the grids we shall derive a 
priori -estimates for and use these to prove that Uh is a quasi-optimal approxima- 
tion of u in the norm in the sense that 

Ik-^i/illi ^ inf \\u-Vh\\i- (6.1) 

While the term on the left in (6.1) is in general difficult to analyze, the term on the 
right represents the fundamental question addressed by classical approximation theory 
in normed spaces, of which much is known. To bound the term on the right from above, 
one picks a function in which is particularly easy to work with, namely, a nodal or 
generalized interpolant of u, and then one employs standard techniques in interpolation 
theory. Therefore, it is clear that the importance of approximation results such as (6.1) 
are that they completely separate the details of the Poisson-Boltzmann equation from the 
approximation theory, making available all known results on finite element interpolation 
of functions in Sobolev spaces (cf. [21]). 

Now we assume can be triangulated exactly (e.g., f2 is a cube) with a shape regu- 
lar and conforming (in the sense of [21]) triangulation Th- Here h = h^^^ represents 
the mesh size which is the maximum diameter of elements in the triangulation. We fur- 
ther assume in the triangulation that the discrete interface Vh approximates the known 
interface V to the second order, i.e., d{T, Th) < Ch'^. 

Given such a triangulation Th of Vt, we construct the linear finite element space V'^ := 
{v e H'^{Q),v\r e Vi{t) \fT e %}. Since the boundary condition g-G e C°^{dQ), 
we can find a uo E H^{^) such that ud = g — G on d^l in the trace sense. Then the 
solution can be uniquely written as u = ud + Uq, with uq E Hq. Thus we will use 
H}j{Q) := Hq^VI) + ud to denote the affine space with a specified boundary condition 
and Vj^ = n Hl){^l) to denote the finite element affine space of H}j{il). Similarly 
Vq = V'^ n i^o (^)- Here to simplify the analysis the boundary condition is assumed to 
be represented exactly. 

Recall that the weak form of RPBE is 

FindM e Hl){n) such that(s.t.) A{u,v) + {B{u),v) + {fG,v) = Vi; G H^{n).(6.2) 

We are interested in the quality of the finite element approximation: 

Find Uh G Vl^ s.t. A{uh, Vh) + {B{uh),Vh) + (/g, v) = Q ^Vh G (6.3) 

It is easy to show that the finite element approximation Uh is the minimizer of E in V^, 
i.e., E{uh) = inf^^gyfe E{vh)- Then the existence and uniqueness follows from section 

3 since is convex. As in the continuous setting, it will be convenient to split the 
discrete solution to the RPBE into linear and nonlinear parts = + u^, where 
and satisfy, respectively. 

Find u\ G s.t. A{ul Vh) + (/g, v) = V^;, G V^, (6.4) 
Find ul G s.t. A«, Vh) + (5« + u'^),Vh) = Vt;^ G V^. (6.5) 
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6.1. Quasi-optimal a priori error estimate. We begin with the following properties of 
the bilinear form A and and operator B. 

Lemma 6.1. 1. The bilinear form A{u, v) satisfies the coercivity and continuity condi- 
tions. That is, for u,v e H^{Q) 

\\u\\\ < A{u,u), and A{u,v) < llMlliUfHi. 

2. The operator B is monotone in the sense that 

{B{u) - B{v),u -v)> r'^Wu - vW^ > 0. 

3. The operator B is bounded in the sense that for u,v & L°°{fl),w e L^(il), 

{B{u) - B{v),w) < C\\u-v\\\\w\\. 

Proof The proof of (1) and (2) is straightforward. We now prove (3). By the mean value 
theorem, there exists 9 e (0,1) such that 

B{u) - B{v) = k'^ cosh{9u + (1 - 9)v + G){u-v). 

Then by the convexity of cosh and the fact that u,v e L°°(Q), G e C°°(Qs), we get 

II cosh(^M + (1 - 9)v + G')||oo,n, < || cosh(M + G')||oo,q, + || cosh(i; + G')||oo,q, < C. 

The desired result then follows since B{-) is nonzero only in fls- D 

Theorem 6.2. Let u and Uh be the solution ofRPBE and its finite element approximation, 
respectively. When Uh is uniformly bounded, we have 

\\u-Uh\\i < inf ||ii-T;/i||i. 

Proof. By the definition, the error u — Uh satisfies 

A{u - Uh, Wh) + {B{u) - B{uh),Wh) =0 Wwhe V^^ 
We then have, for any Vh e V^, 

WtJ'-tJ'hWl ^ A{u-Uh,u-Uh)^A{u-Uh,u-Vh) + A{u-Uh,Vh-Uh) 
< \\u - Uh\\i\\u - Vh\\i - iB{u) - B{uh),Vh - Uh). 

The second term on the right side is estimated by 

-{B{u) - B{uh),Vh-Uh) = -iB{u)- B{uh),u-Uh) + {B{u)-B{uh),u-Vh) 

< {B{u)-B{uh),u-Vh) 

^ \\u - Uh\\l\\u - VhWl- 

Here we make use of the monotonicity of B in the second step and the boundness of B 
in the third step. In summary we obtain for any Vh £ 

Ik - Uh\\l < \\U - Vh\\l, 

which leads to the desired result by taking the infimum. □ 
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Figure 1 . Divide a cube into 5 tetrahedra. 

6.2. Discrete a priori -estimates. We now derive L°°-estimates of the finite element 
approximation Uh. To this end, we have to put assumptions on the grid. Let (aij) denote 
the matrix of the elliptic operator (eVw, Vv), i.e., = A{ipi, ipj). Two nodes i and j 
are adjacent if there is an edge connecting them. 

(Al) The off-diagonal term aij, i,j are adjacent, satisfies 

We now give example grids satisfying (Al). In three dimensions, to simplify the 
generation of the grid, we choose as a cube and divide into small cubes with length 
h. For each small cube, we divide it into 5 tetrahedra; see Figure 1 for a prototype of 
the triangulation of one cube. Neighbor cubes are triangulated in the same fashion (with 
different reflection to make the triangulation conforming). By the formula of the local 
stiffness matrix in [32, 54], it is easy to verify that the grids will satisfy assumption (Al). 
We comment that the uniform grid obtained by dividing each cube into 6 tetrahedra will 
not satisfy the assumption (Al), since in this case if i, j are vertices of diagonal of some 
cube, then aij = 0. 

Theorem 6.3. In general dimension W^,d >2, with assumption (Al) and h sufficiently 
small, the finite element approximation Uh ofRPBE satisfies 

\\Uh\\oo ^ C, 

where C is independent ofh. 

Proof. We shall use the decomposition u/i = + u^. By the regularity result [41], we 
know G bI^^{Q,) and thus obtain a priori estimate on quasi-uniform grids 

II'"' - "illoo < Ch'^^^ < Cdiam(fi)^ for some s E (0, 3/2). 

This implies that Hm^Hoo < ||w'||oo + ||^^' — "J^/iHoo < C is uniformly bounded with respect 
to haiax- The estimate of follows from Theorem 3.3 in [32], where the grid assumption 
(Al) is used. □ 

In two dimensions, we can relax the assumption on the grid and obtain a similar result. 
Later we will see that, due to this relaxation, the local refinement in two dimensions is 
pretty simple. 

(Al') The off-diagonal terms j < 0, j 7^ i; i.e., the stiffness matrix corresponding to 
■) is an M-matrix. 

Theorem 6.4. For a two-dimensional triangulation satisfying (Al'), the finite element 
approximation Uh ofRPBE is bounded, i.e., 

\\Uh\\oo < C. 
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Proof. Similarly ||m/j||oo < C is uniformly bounded. In two dimensions the estimate of 
n^l follows from Theorem 3.1 in [32], where the grid assumption (Al) is 
used. □ 



7. Convergence of adaptive finite element approximation 

In this section, we shall follow the framework presented in [49, 50] to derive an a 
posteriori error estimate. Furthermore we shall present an adaptive method through local 
refinement based on this error estimator and prove that it will converge. The a priori 
L°° -estimates of the continuous and discrete problems derived in the previous sections 
play an important role here. 

7.1. A posteriori error estimate. There are several approaches to adaptive error con- 
trol, among which the one based on a posteriori error estimation is usually the most ef- 
fective and most general. Although most existing work on a posteriori estimates has been 
for linear problems, extensions to the nonlinear case can be made through linearization. 
For example, consider the nonlinear problem 

F(n) = 0, FeC\B^,Bl), (7.1) 

where the Banach spaces Ei and B2 are, e.g., Sobolev spaces and where B* denotes the 
dual space of B. Consider now also a discretization of (7.1) 

Fh{uh) = ^, FneC\Uh,V:), {12) 

where Uh C Bi and Vh C B2- For the RPBE and a finite element discretization, the 
function spaces would be taken to be Bi = B2 = HliVt). The nonlinear residual F{uh) 
can be used to estimate the error through the use of a linearization inequality 

Ci||FK)||b- < \\u-Uh\\B, < C2\\F{uh)\\B',. (7.3) 

See, for example, [49] for a proof of this linearization result under weak assumptions on 
F. The estimator is then based on an upper bound on the dual norm of the nonlinear 
residual on the right in (7.3). 

In this section, to show the main idea, we will assume Fh{uh) = F{uh) by making the 
following assumption on the grid. 

(A2) The smooth interface V is replaced by its discrete approximation Vh such that e 
and K are piecewise constants on each element of the triangulation 7/, . 

In our setting of the weak formulation, we need to estimate \\F{uh) To this end, 
we first introduce quite a bit of notation. We assume that the rf-dimensional domain Vt has 
been exactly triangulated with a set Th of shape -regular d-simplices (the finite dimension 
d is arbitrary, not restricted to c? < 3, throughout this discussion). A family of simplices 
will be referred to here as shape-regular in the sense of [21]. 

It will be convenient to introduce the following notation: 





= the set of shape-regular simplices triangulating the domain Q. 


M{r) 


= the union of faces contained in simplex set r lying on dfl. 


X(r) 


= the union of faces contained in simplex set r not in N'ir). 


Hr) 


= U{t)\JX{t). 


J" 


= UrenHr)- 


UJr 


= [j {f eTh\T C\f ^ 0, where T eTh}. 


Us 


= [J{ferh\Sf]fy^0, where SeT}. 




= the diameter of the simplex r. 


hs 


= the diameter of the face S. 
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When the argument to one of the face set functions A/", X, or T is in fact the entire set of 
simplices, we will leave off the explicit dependence on S without danger of confusion. 
Finally, we will also need some notation to represent discontinuous jumps in function 
values across faces interior to the triangulation. For any face S E M, let ns denote the 
unit outward normal; for any face S E X, take ns to be an arbitrary (but fixed) choice 
of one of the two possible face normal orientations. Now, for any v E L^(fi) such that 
V E C°(r) Vr E Th, define the jump function: 

\v]s(x) = lim v(x + tns) — lim vix — tns). 
We now define our a posteriori error estimator 

ril{uh) = hl\\B{uh)\\l, + ]^ J2 hs\\[ns ■ {eVuh + {e - ejVG)]s\\ls, (7.4) 

5GX(r) 

and the oscillation 

osc^(m,) = K {W'^UhWlr + llVGIIt) . (7.5) 

Theorem 7.1. Let u E H^{Q) be a weak solution of the RPBE and Uh be the finite 
element approximation with a grid satisfying assumptions (Al) and (A2). There exist 
two constants depending only on the shape regularity of Th such that 

\\u-Uh\\\<Cirjl + C20scl, (7.6) 

where 

Proof. We shall apply the general estimate in [50, Chapter 2] (see also [49]) to 

a{x,u,'Vu) = eVu + {e — EmJ'VG, and b{x,u,'Vu) = —r'^ smh.{u + G). 

We then use the following facts to get the desired result: 

• V ■ {eVuh) |r = Vr G by the assumption (A2) of the grid; 

• V ■ ((e - £„)VG) |r = Vr G Tfc since AG(x) = if x ^ {xi}. 

• For r G Th flf^s, let Uh and G denote the average of Uh and G over r, respectively. 
We then have 

II sinh(M/i + G) - sinh(M/, + G)||o,r < \cosh.{0\\\uh - Uh + G - GW^^r 

< C/i2(||Vm,||o,. + ||VG||o,.). 

Here we use the L°°-estimates of u and Uh to conclude that | cosh(^)| < G and 
the standard error estimate for Hm/j — -u/iUcr and ||G — (5||o,r- 

We give some remarks on our error estimator and the oscillation term. First, using 
(4.2) one can easily show that ||VM/i||o,f7 < G uniformly with respect to h and thus 
osCr = 0{hl). Comparing to the order of r]r = 0{hr), the error estimator r/^ will 
dominate in the upper bound. Second, in (7.4) the jump of [ns ■ {e — VG]^ 7^ only 
as E Vh- This additional term with order 0([£:]) will emphasize the elements around 
the interface where the refinement most occurs. 

Although it is clear that the upper bound is the key to bounding the error, the lower 
bound can also be quite useful; it can help to ensure that the adaptive procedure does 
not do too much work by overrefining an area where it is unnecessary. Again using the 
general framework for the a posteriori error estimate in [49, 50], we have the following 
lower bound result. 
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Theorem 7.2. There exists two constants C3, C4 depending only on the shape regularity 
ofTh such that 

vliuh) < CsWu - Uh\\l^^^ + C4 ^ oscl{uh) Vr G %■ 

7.2. Marking and refinement strategy. Given an initial triangulation To, we shall gen- 
erate a sequence of nested conforming triangulations Tk using the following loop: 

SOLVE ESTIMATE MARK REFINE. (7.7) 

More precisely to get T^+i from Tk we first solve the discrete equation to get Uk on 
Tfc. The error is estimated using Uk and used to mark a set of triangles that are to be 
refined. Triangles are refined in such a way that the triangulation is still shape-regular 
and conforming. 

We have discussed the step ESTIMATE in detail, and we shall not discuss the step 
SOLVE, which deserves a separate investigation. We assume that the solutions of the 
finite-dimensional problems can be solved to any accuracy efficiently. Examples of such 
optimal solvers are the multigrid method or the multigrid-based preconditioned conjugate 
gradient method [52, 13, 27, 55]. In particular we refer to [1, 2] for recent work on 
adaptive grids in three dimensions and [30, 29] for solving the PBE with inexact Newton 
methods. 

We now present the marking strategy which is crucial for our adaptive methods. We 
shall focus on one iteration of loop (7.7) and thus use Th for the coarse mesh and Th for 
the refined mesh. Quantities related to those meshes will be distinguished by a subscript 
H or h, respectively. 

Let 9i, i = 1, 2 be two numbers in (0, 1). 

(1) Mark Mi^h such that 

reMi,H reTn 

(2) If 

OSCh > TjH (7.8) 

or 

C, 0SC2(W^) >\Y. "^rM^ (7-9) 

then extend M.i^h to M.2,h such that 

reM2,H reTn 

Unlike the marking strategy for reducing oscillation in the adaptive finite element meth- 
ods in [36, 37], in the second step, we put a switch (7.8)-(7.9). In our setting, the oscil- 
lation OSCh = 0{H'^) is in general a high-order term. The marking step (2) is seldom 
applied. 

In the REFINE step, we need to carefully choose the rule for dividing the marked 
triangles such that the mesh obtained by this dividing rule is still conforming and shape- 
regular. Such refinement rules include red and green refinement [7], longest refinement 
[40, 39], and newest vertex bisection [42, 34, 35]. For the REFINE step, we are going to 
impose the following assumptions. 

(A3) Each r G Aiu, as well as each of its faces, contains a node of Th in its interior. 
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(A4) Let 7h be a refinement of Th such that the corresponding finite element spaces 
are nested, i.e., C V^. 

With those assumptions, we can have the discrete lower bound between two nested 
grids. Let 7}/ be a shape-regular triangulation, and let Th be a refinement of Th obtained 
by local refinement of marked elements set M.h- The assumption (A3) is known as the 
interior nodes property in [37]. Such a requirement ensures that the refined finite element 
space V'^ is fine enough to capture the difference of solutions. 

Theorem 7.3. Let Th be a shape-regular triangulation, and let Th be a refinement of 
Th obtained by some local refinement methods of marked elements set Ai h, such that it 
satisfies assumptions (A3) and (A4). Then there exist two constants, depending only on 
the shape regularity ofTn, such that 



Proof. The proof is standard using the discrete bubble functions on r and each face 



7.3. Convergence analysis. We shall prove that the repeating of loop (7.7) will produce 
a convergent solution Uk to u. The convergent analysis of the adaptive finite element 
method is an active topic. In the literature it is mainly restricted to the linear equa- 
tions [17, 46, 16, 36, 25, 9, 37, 33, 26, 8]. The convergence analysis for the nonlinear 
equation is relatively rare [24, 48, 15]. 

Lemma 7.4. Let Th and Th satisfy assumptions (A1)-(A4). Then there exist two con- 
stants depending only on the shape regularity ofTn such that 



rjliun) < C^\\uh - Uh\\\u,^ + '^4 ^ osc?(m/^) Vr G Mh- 



(7.10) 



S e dr. 



□ 



\\u - uhWI < C^\\uh - uhTi + C'eosc^. 
When (7.8) and (7.9) do not hold, we have a stronger inequality 

\u - unfi < Cr \\uh - uh\\1, 



where C-j depends only on the shape regularity ofTn- 
Proof. By the upper bound and marking strategy 



u-uh\\\ < Cir]jj + C20sc% 



T&Ml,H 




Cs = Ci^r C's" , and Cg = (C2 + 2C^^Ci). 
If (7.8) does not hold, i.e., osch < ?7h, the first inequality becomes 

\\u-uh\\1 < {Ci + C2)Vh- 
If (7.9) does not hold, we can easily modify the lower bound (7.10) as 




Then the inequality follows similarly. 



□ 
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For T/i C th, let /i^^^ = 'jHt-^, with 7 G (0, 1). The next lemma shows that even the 
oscillation is not small; there is also a reduction result. For the marked set A^// C Th, 
we shall use to denoted the refined elements in Th- 

Lemma 7.5. If M.2,h\M.i,h ^ 0, there exist pi, P2 such that 

OScl < Pi OScjj + p2\\Uh - Uh\\1- 



Proof. 

< oscI{uh)+ Y1 oscl{uH) + Ch^\\V{uh-UHW 

< 7^ ^ oscI{uh)+ oscl{uH) + Ch''\\V{uh-UHW 

TH&M2,H THeTH\M2,H 

< os4 + (7'-l) osclM + Ch^Viuh-unW 

TH&M2,H 

< PlO?,c]j + p2\\Uh- Uh\\\, 

with pi = 1 - (1 - 72)/^2 e (0, 1), and p2 = CK^. □ 

We shall choose 62 sufficiently close to 1 and /imax < 1/c to ensure pi G (0, 1), i = 
1,2. 

For the nonlinear problem, we do not have the orthogonality in norms. But we 
shall use the trivial identity 

E{uh) - E{u) = E{uh) - E{uh) + E{uh) - E{u). (7.11) 

The following lemma proves the equivalence of energy error and error in norm. Again 
the norm estimate of u and Uh is crucial. 

Lemma 7.6. If both Th and Th satisfy the assumption (Al), then 
• E{uh) - E{u) ~ \\uh - u\\l; 



E{uh) — E{u) ~ \\uh — u 



|2. 

1' 



• E{uh) - E{uh) ^ \\uh - UhWl- 

Proof. By the Taylor expansion 

E{uh) - E{uh) = {DE{uh),UH - Uh) + {D'^E{^){uh - Uh), uh - Uh). 

The first term is zero since Uh is the minimizer. The desired result follows from the bound 

< \\D''E{0\\oo = R^W cosh(C + G')||oo,n. < C. 

Other inequalities follow from the same line. □ 

Our adaptive finite element methods (AFEMs) consist of the iteration of loop (7.7) 
with the estimate, marking, and refinement parts discussed before. Also the grids gener- 
ated by the algorithm will satisfy assumptions (A1)-(A4). Hereafter we replace the sub- 
script h by an iteration counter called k and introduce some notation to simplify the proof. 
Let Uk be the solution in the kth iteration, 5k := E{uk) — E{u), dk = E{uk) — E{uk+i), 
and Ok = osc^('Ufc) 
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Theorem 7.7. The adaptive method using loop (1.1) will produce a convergent approxi- 
mation in the sense that 

lim \\u — UkWi = 0. 

Proof. By Lemma 7.6, we need only to show 6'' as k ^ 0. We first discuss tlie 
easier case: When gsc^ is the high-order term in the sense that the inequalities (7.8) and 
(7.9) do not hold, we have the error reduction 

\\u - uhWI < C\\uh - uhW^. 

Using Lemma 7.5 and (7.1 1), we have 

E{uh) - E{u) < C{E{uh) - E{uh)), 

which is equivalent to 5h < C5h — C5h. Then 6h < (1 — 1/C)5h, and thus 

6'' < a^5^, with a = (1 - l/C) e (0, 1). 

When the oscillation is not small, i.e., (7.8) or (7.9) holds, we can get only 

Ai4 <dk + A20A,., with Ai G (0, 1). (7.12) 

We shall use techniques from [33] to prove the convergence. Recall that we have 

4+1 = Sk - 4- (7.13) 

For any /3 e (0, 1), /3 x (7.12) + (7.13) gives 

5fc+i < «4 + /3A20fc - (1 - /3)4, witha = (l-/3Ai) G (0,1). (7.14) 

Recall that we have 

Ofe+i < PiOk + P2dk. (7.15) 
Let 7 = (1 - /5)/p2; (7.15) X 7 + (7.14) gives 

Sk+i + lOk+i < ah + (/3A2 + pi7)ofc- 
Let 1 > /i > pi . We choose 

to get 

5k+i + 70fc+i < max(a, ^){5k + 70fc), 
which also implies the convergence of our AFEM. □ 

8. Summary and concluding remarks 

In this article we have established a number of basic theoretical results for the nonlin- 
ear Poisson-Boltzmann equation and for its approximation using finite element methods. 
We began by showing that the problem is well-posed through the use of an auxiliary or 
regularized version of the equation and then established a number of basic estimates for 
the solution to the regularized problem. The Poisson-Boltzmann equation does not ap- 
pear to have been previously studied in detail theoretically, and it is hoped that this paper 
will help provide molecular modelers with a better theoretical foundation for their an- 
alytical and computational work with the Poisson-Boltzmann equation. The bulk of 
this article then focused on designing a numerical discretization procedure based on 
the regularized problem and on establishing rigorously that the discretization procedure 
converged to the solution to the original (nonregularized) nonlinear Poisson-Boltzmann 
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equation. Based on these results, we also designed an adaptive finite element approxi- 
mation procedure and then gave a fairly involved technical argument showing that this 
adaptive procedure also converges in the limit of mesh refinement. This article appar- 
ently gives the first convergence result for a numerical discretization technique for the 
nonlinear Poisson-Boltzmann equation with delta distribution sources, and it also intro- 
duces the first provably convergent adaptive method for the equation. This last result is 
one of only a handful of convergence results of this type for nonlinear elliptic equations 
(the others being [24, 48, 15]). 

Several of the theoretical results in the paper rest on some basic assumptions on the 
underlying simplex mesh partitioning of the domain, namely, assumptions (A1)-(A4); 
we now make a few comments on these assumptions. To begin, we required a refinement 
procedure that would preserve the norm estimate of Uh. Meeting this requirement 
in the two-dimensional setting is relatively easy; one can choose as a square and start 
with a uniform mesh of a square. For the refinement methods, one can use longest edge 
or newest vertex bisection. Subdivisions obtained by these two methods contain only one 
type of triangle: isosceles right triangles. Thus the assumption (Al') always holds. In 
the three-dimensional setting, this is more tricky. Bisection will introduce some obtuse 
angles in the refined elements. One needs to use a three-dimensional analogue of red- 
green refinement [10]. However, this will not produce nested subspaces; i.e., assumption 
(A4) is invalid. For convergence analysis based on red-green refinement, we could use 
the technique in [45] to relax the assumption (A4). Since this will only add technical 
difficulties but does not exhibit principally new phenomena, we omit them here. Another 
approach to relax the assumption (Al) is to use pointwise a posteriori error estimates 
developed in [38] for monotone semilinear equations. We can start with a quasi-uniform 
triangulation and refine the triangulation according to the pointwise a posteriori error es- 
timator to make sure \\u — Uh\\oo < C. Then together with the L°° norm estimate of u, 
by the triangulation inequality ||Mh||oo < ||m||oo + ||w — m/i||oo < C, we have the control 
of ||m/i||oo- Note that the pointwise a posteriori error estimates developed in [38] are for 
elliptic-type equations with continuous coefficients. To use this approach we need to 
adapt the estimate for the jump coefficients case which will be a further research topic. 

Assumption (A2) is needed to approximate the interface well in an a priori manner. 
Of course, one can include this approximation effect into the a posteriori error estimate 
(namely, the term \\F(uh) — Fh{uh)\\) and use this to drive local refinement to improve 
the approximation to the desired level for the assumption or use the strategy for the os- 
cillation to include it in the refinement loop. However, we note that, since the interface is 
known a priori from, e.g., x-ray crystallography information, we do not need to solve the 
equation (which is generally the more expensive route) to solve this problem; we view 
this as primarily a mesh generation problem. Robust algorithms to produce well-shaped 
tetrahedral meshes which are constrained to exactly match some interior embedded two- 
manifold are available in the literature; for example, see [18, 3]. A simple algorithm 
can be based entirely on local refinement with the marking and refinement strategy, but 
without having to solve the PBE to produce error indicators: If the element cross the 
interface, then it gets refined. This strategy was employed in [5]. 

After this work was done, we learned that the assumption (A3) is not needed for the 
convergence of adaptive finite element methods for a linear elliptic equation. As an 
ongoing project, we are extending it to the nonlinear Poisson-Boltzmann equation. 

Finally, we make some remarks on the practical realization of a convergent discretiza- 
tion procedure based on the two-way (or three-way) expansion into a known singular 
function and solution(s) of an associated regularized version of the problem. Methods 
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for building high-quality approximate solutions of the regularized nonlinear PBE, either 
by solving (3.5)-(3.6) at once or by solving for the linear and nonlinear pieces separately 
by solving (5.1)-(5.2) and then adding the solutions together, are well-understood. The 
techniques described in [28], taken together with the approximation framework and the 
adaptive algorithm proposed in the present article, moves us a step closer to the goal 
of a complete optimal solution to this problem, in terms of approximation quality for 
a given number of degrees of freedom, computational complexity of solving the corre- 
sponding discrete equations, and the storage requirements of the resulting algorithms. 
What remains is simply the cost of evaluating the singular function G in forming the 
source terms in (3.5) or (5.1). The source terms are evaluated using numerical quadra- 
ture schemes: sampling the integrand at specially chosen discrete points in each element 
and then summing the results up using an appropriate weighting. This is equivalent to 
computing all pairwise interactions between the collection of quadrature points (a fixed 
constant number of points per simplex) and the number of charges forming G. Given that 
G is typically formed from at most a few thousand charges, the algorithm evaluating G at 
the quadrature points should scale linearly with the number of quadrature points, which 
is a (small) constant multiple of the number of simplices. This can be accomplished 
using techniques such distance-classing and fast multiple-type methods. 
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